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Study of turbulent premixed flame 
propagation using a laminar flamelet model 

By H. G. Im 


1. Motivation and objectives 

The laminar flamelet concept in turbulent reacting flows is considered applicable 
to many practical combustion systems (Linan h Williams 1993). For turbulent 
premixed combustion, the laminar flamelet regime is valid when turbulent Karlovitz 
number is less than unity, which is equivalent to stating that the characteristic 
thickness of the flame is less than that of a Kolmogorov eddy; this is known as the 
Klimov-Williams criterion (Williams 1985). In such a case, the flame maintains its 
laminar structure, and the effect of turbulent flow is merely to wrinkle and strain 
the flame front. The propagating wrinkled premixed flame can then be described 
as an infinitesimally thin surface dividing the unburnt fresh mixture and the burnt 
product. 

It has been suggested (Kerstein et al. 1988) that such a propagating front can be 
represented as a level contour of a continuous function G, whose governing equation, 
derived using the Huygens’ principle, is 

'§ + «s;-« |vc i- (1) 

Here sl is the well-defined laminar flame speed which is generally not a constant, 
but can be modified by the effect of flame stretch. By introducing the Markstein 
length C (Pelce & Clavin 1982), an asymptotic analysis gives an expression for si\ 

s L - s° L - s° L CV n + Cn • (Vn) * n, (2) 

where n = — VG/|VG| is the normal vector to the surface pointing toward the 
unburnt mixture. The Markstein length is of the order of flame thickness A / pCpSi 
defined usually in terms of unburnt mixture properties. Here A is the thermal 
conductivity and c p the specific heat. 

There are several advantages to using the G-equation model rather them direct 
numerical simulation with Arrhenius-type chemistry. First, since the flame front 
is described by a contour of the smooth function G, complex topology changes 
in the propagating front can be easily captured by solving the transport equa- 
tion for G, instead of tracking the corrugated front. Secondly, since the numerical 
stiffness due to the Arrhenius chemistry with large activation energy is removed 
in favor of a flamelet whose structure is given a priori , the computational effort 
can be significantly reduced with an appropriate discontinuity-capturing numerical 
scheme. Furthermore, the diffusional-thermal modification of the flame structure is 
accounted for by the flame-speed relation (2) in a parametric manner; the coupling 
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between the hydrodynamic field and the flame-structure is simply accounted for by 
the parameter C . This is important in validating the existing predictions of turbu- 
lent flame speed, most of which are based on the constant sl assumption. Finally, 
by eliminating the nonlinear reaction terms from the conservation equations, the 
system can be more easily adapted to large-eddy simulation based on the dynamic 
subgrid-scale modeling principle. A preliminary attempt at such modeling will be 
discussed in a later section. 

From a fundamental standpoint, the G-equation model serves as a useful tool for 
understanding some issues in turbulent premixed combustion. One such issue is the 
determination of turbulent flame speed, sr, as a function of flow quantities such as 
the turbulence intensity, v! . Although there are theoretical models and experimental 
observations, the agreement among the various results is far from being satisfactory. 
Thus far, perhaps the only concensus is that st increases with u f initially, then tends 
to level off at larger id, which is often called “bending” behavior (Bradley 1992). 

Most theoretical models of st in the flamelet regime are based on Damkohler’s 
(1940) proposition that the increase in the flame speed is proportional to the area 
increase, which in turn can be related to the turbulence intensity. This suggests 

st/$l = At/Al ~ 1 + C(u' fsi) q , (3) 

where At is the total surface area of the wrinkled front and Al the cross-section 
area normal to the direction of propagation. Based on this proposition, Clavin 
& Williams (1979) derived q = 2 from geometrical considerations, while Yakhot’s 
renormalization group theory (1988a) yields the same result in the weak turbu- 
lence limit. Recently, Kerstein &: Ashurst (1992) proposed q = 4/3 by considering 
the random nature of turbulent flows. This result was further supported by their 
numerical study (Kerstein & Ashurst 1994). 

All of these arguments are based on the constant density assumption so the effect 
of heat release generated by chemical reaction has not been taken into account. 
Variable density introduces additional complexities, one being that the coupling 
between flow and flame must be dealt with. Recently, Cambrav & Joulin (1992), 
in a semi-analytic study of the model equation by Michelson & Sivashinsky (1977), 
demonstrated that, at least if u f < 0 (sl ), the turbulent burning velocity is no- 
ticeably enhanced by hydrodynamic instability. Their numerical results suggest the 
value q of about 0.3 in the weak turbulence range. If validated by further studies, 
this result may show that the “bending” behavior may be the effect of thermal- 
expansion induced wrinkling, which diminishes at higher u'. 

Therefore, in this study we attempt to provide a useful database for understanding 
these issues in turbulent premixed combustion. In particular, the effect of thermal 
expansion is investigated by fully coupling the G-equation with the flow field. In the 
following section, the formulation of the variable- density version of the G-equation 
model is presented, and some numerical results are discussed for premixed flames 
propagating in a harmonic inlet velocity flow field and a pair of counter-rotating 
vortices. The results of the former problem are consistent with those of Cambray &; 
Joulin (1992), while the study of the flame- vortex interaction also reveals interesting 
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behaviors regarding the vorticity produced by flame. Finally, a subgrid-scale model 
for the G-equation based on the dynamic modeling concept is proposed. 

2. Accomplishments 

2.1 The G-equation model with heat release 


2.1.1 Formulation 

Throughout this study, we define the flame front as the contour, G = 0, of a 
continuous function G(x,t), where G < (>) 0 is defined as the unburnt (burnt) 
side. The species equation is then substituted by the G-equation which can be 
written in conservative form as (Williams 1985) 


I ( ' ,G) + S7 < ' ra ' G, = '“ t|vc| - 


( 4 ) 


Using the flame-speed relation (2) with the definition n = — V(r/|VG|, we obtain 
(Peters 1992) 


I <*0 ■ + KJ <««> - ** (| + GV ’ G - If) 


1 duk dG dG 
+pC \ VG\dxjdx]d^ 


( 5 ) 


where the subscript 0 denotes the condition at the unbumt mixture, s° L the plane 
laminar flame speed, and we use the approximation psi = PqS° l — constant. Equa- 
tion (5) accounts for the effect of the flame stretch given by the results (2). 

To include the effect of thermal expansion, we introduce the total energy 


e = l u 2 i+ c v T + q[l-H(G)] (6) 

where H is the Heaviside function. This implies that as the flow crosses the flame 
(G = 0), an amount of chemical energy q is converted to thermal energy, thereby 
creating jumps in the density and temperature. The conservation equation for the 
total energy is free of reaction term, i.e. 


m {pc) + (0,e + p, “' 1 



( 7 ) 


where P is the pressure, r,j the stress tensor, and the heat flux is given by Fourier’s 
law. 

The rest of the system consists of the continuity equation 


dp d , x „ 


( 8 ) 
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the momentum equation 


'} 


d , . d dP dr, 


( 9 ) 


and the equation of state 


P = pRT. 


( 10 ) 


In the present numerical simulations, the discontinuity is removed by replacing the 
Heaviside function by the smooth function 


H{G) « [1 + tanh(G/6|VG|)], (11) 


where 8 is a small parameter of the order of the flame thickness. 

The fully-compressible system (5)-(10) is solved using a high order compact 
scheme (Lele 1992) for spatial derivatives and a third order Runge-Kutta scheme 
(Wray 1990). Boundary conditions axe treated following the method of Poinsot 
and Lele (1992). For one-dimensional calculations, the initial condition for the G 


function is 


G(x) 


' -1, if x — Xf < —W\ 

i sin[n(x — x f)/2W], if |x — x/| < W; 
1, if x — Xf > W 7 


( 12 ) 


and the boundary condition on G is treated in the same way as the other scalar 
variables. Here W is the thickness of the G profile. The converged one- dimensional 
solution is used as the initial condition for the two-dimensional calculation. 

Figure 1 shows schematics of the two model problems considered, namely the 
flame response to (a) a steady harmonic velocity fluctuation, and (b) a pair of 
counter-rotating vortices. Some results for each model problem are presented and 
discussed below. 


2.1.2 Harmonic inlet velocity 

As shown in Fig. 1(a), we impose a steady harmonic inlet velocity profile 


u(x = 0, t) = s° L + u cos(27ry). 


(13) 


For u l = 0, the G-field remains fixed at the initial condition. In a simulation, at 
t = 0 a finite value of u f is imposed at the inlet boundary; this velocity fluctuation 
then produces a curved front. The calculation proceeds until a final state is attained, 
in which the flame area does not change and the front moves toward the unburnt 
mixture due to the enhanced propagation rate. In the present calculation we used 
the parameter values Re a = [aLj i/) 0 = 2000, where a is the speed of sound, unity 
for the Prandtl and Lewis numbers, and s° L /a = 0.05. The results depend on the 
Markstein length C through the flame-speed variation (see (2)). To minimize this 
flame-structure effect and to extract the behavior of the flame in the Huygens’ limit, 
we choose CjL = 0.01 in the present calculation, where L is the width of the channel 
shown in Fig. 1(a). 
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u(x=0) = u 0 + cos(2xcy/L) Periodic 
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Nonreflecting 
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FIGURE 1. Schematics of the model problems and computational conditions: (a) 
flame propagating into the steady harmonic inlet velocity, (b) flame-vortex inter- 
action, where the solid and dotted curves respectively denote positive and negative 
vorticities. 


Figure 2 shows the final state of the flame fronts represented by the G = 0 
contours for the inlet perturbations of u f /s° L = 0 and 0.3. Here a = ( p u — Pb)/Pu is 
the heat release parameter; a = 0 for the zero heat-release case and a = 0.5 when 
the downstream temperature is twice the upstream temperature. It is seen that 
the flames with heat release (a = 0.5) are more curved than those without heat 
release (a = 0). This is due to the hydrodynamic instability mechanism known as 
the Landau-Darrieus effect (Williams 1985). At a hydrodynamic discontinuity with 
constant propagation speed, thermal expansion induces a deflection of streamlines 
such that the convex front is further accelerated. Although the linear stability 
analysis predicts that the perturbation of the front grows indefinitely, in reality 
it saturates as nonlinear effects come into play. Figure 2 clearly demonstrates 



352 


H. G. Im 



FIGURE 2. Flame fronts described as G = 0 contours subject to the steady 
harmonic inlet velocity for (a) u 1 /s° L = 0 and (b) u'/s° L = 0.3. Shown in each figure 
are the cases for zero heat release (a = 0) and for a = 0.5. 

such behavior, and the flame propagating with larger heat release exhibits more 
wrinkling. In particular, it is of interest to note from Fig. 2(a) that with heat 
release the flame front does not remain planar even if inlet velocity perturbation is 
absent ( u r — 0), consistent with the result of Cambray & Joulin (1992). 

In Fig. 3 we plot the area ratio ( At/Al ) as a function of the magnitude of 
velocity fluctuation (v* /$° L ) for the configuration shown in Fig. 1(a). At present, the 
range of u* /s° L is limited due to numerical difficulty that arises when u ' significantly 
exceeds S° L SO that the front forms sharp curvature. Nevertheless, Fig. 3 confirms 
the results of Cambray & Joulin (1992) in that there is an additional flame-speed 
enhancement due to thermal expansion for weak turbulence (u' /s° L < 1). For larger 
velocity fluctuations, it is expected that the effect of thermal expansion induced 
self-wrinkling of the front will be less prominent as the large convective flow field 
dominates the flame behavior, which may be a possible mechanism for the “bending” 
behavior. Further improvement in the numerical methodology to capture more 
excessive wrinkled front is required to obtain a more conclusive database regarding 
this issue. 

2.1.3 Flame-vortex interaction 

To further investigate the coupling between a flame and a flow via density varia- 
tion, we adopt the flame-vortex interaction as a model problem, as wa s previously 
studied by Poinsot et ai (1991). In particular, the emphasis is on fundamental 
issues such as the flame front response to the vortical flows and attenuation and 
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FIGURE 3. Nondimensional total flame-front area vs. nondimensional velocity 
fluctuation. Open symbols are for zero heat release (a = 0) and solid symbols for 
a = 0.5. 

generation of vorticity by the flame due to thermal expansion. As shown in Fig. 1(b), 
at t = 0 we introduce a pair of counter-rotating vortices into the uniform flow field 
with uo = s° L , far upstream of the flame. Then, due to the mean flow as well as the 
flow induced by the vortices, the vortex pair drifts downstream and passes through 
the propagating flame front, while preserving symmetry. The initial circulation, I\ 
of the vortices adopted in this study is given by 

r(r) = ±2tt^^2 exp (“^ 2 ) » ( 14 ) 

where r is the distance from the vortex center and a the characteristic radius of the 
vortex. Here we define the strength of the vortex v! by the maximum circumferential 
velocity at t = 0. Other parameter values used in this study are Re a = 1000, 
Pr = Sc — 0.75, fi/fio = (T/To) 0 - 76 , s° L /a = 0.02, C/(X/pc p s° L ) 0 = 0.1. The vortex 
diameter is initially about three times larger than the flame thickness and grows in 
time by diffusive transport. 

Figures 4 and 5 show the snapshots of the flame front and vorticity contours at 
the instant that the flame is most wrinkled by the vortex, for two vortex strengths, 
(u'/s ° L ) t = 0 = 2.4 and 4.8. In each figure, (a) is for the cold flame case (a = 0) 





354 


H. G. Im 



(a) (b) 

FIGURE 4. Flame-vortex interaction for u' /s° L = 2.4, (a) a = 0 and (b) a = 0.75. 
Top and bottom figures respectively denote flame fronts (G = 0) and vorticity. The 
solid and dotted curves respectively denote counterclockwise and clockwise vorticity. 


and (b) for a = 0.75. Although not presented here, the results of the G-equation 
model have been compared to that with the one-step Arrhenius chemistry, and it 
was found that the G-equation captures the essential physics of the flame and flow 
responses. It is also remarked that, due to the rapid decrease in the tangential 
velocity for the initial field (14), an additional vortex pair with opposite sign is 
formed behind the incident vortex pair. Although it may be unphysical, this fast- 
decaying vortex requires a smaller computational domain, and thus adopted in this 
qualitative study. 

Figure 4 is for the lower vortex strength. It is seen that, while the vortices 
Fig. 4(a) preserve their original shapes through the flame, in Fig. 4(b) the vortices 
are significantly elongated behind the flame due to thermal expansion accelerating 
the flow. Furthermore, in this case it is interesting to note that the sign of the 
vorticity is reversed as the vortex passes through the flame. This demonstrates the 
vorticity generation due to the baroclinic torque mechanism arising from the fact 
that the pressure and density gradients are not parallel across the curved flame. In 
this configuration the flame-generated vorticity is opposite to the incident vorticity. 
Therefore, for the case shown in Fig. 4(b), the incident vortices is overridden by the 
flame- generated vortices and cannot survive the flame. Consequently, the reversed 
velocity field induced by the flame-generated vorticity tends to push the retarded 
flame front forward, yielding a less wrinkled front compared to the cold-flame case 
shown in Fig. 4(a). The results agrees qualitatively with a recent experimental 
observation (Mueller et al. 1995). 

Figure 5 shows the case of a stronger vortex, u* /s° L = 4.8. The front becomes 
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(a) (b) 


FIGURE 5. Flame- vortex interaction for u* /s° L = 4.8, (a) a = 0 and (b) a = 0.75. 
Top and bottom figures respectively denote flame fronts (G = 0) and vorticity. The 
solid and dotted curves respectively denote counterclockwise and clockwise vorticity. 

more wrinkled. Consistent with the results in Fig. 4, it is seen that the flame-front 
wrinkling is less severe in the a = 0.75 case. Unlike Fig. 4(b), however, the inci- 
dent vortices shown in Fig. 5(b) axe sufficiently strong to survive the flame, except 
around the sharply curved front where the flame-generated vorticity is most intense. 
Although the vorticity downstream of the flame has the same sign as the incident 
vorticity, the strength of the vorticity is considerably weakened. The mechanisms 
of the vorticity attenuation by the flame are the aforementioned flame-generated 
vorticity and volume expansion, which spreads out the vortical region while pre- 
serving the total circulation (cf. Mueller et al. 1995). These front-stabilizing effects 
may be partly responsible for the experimentally observed “bending” behavior of 
st at high turbulence levels, along with the hydrodynamic effect discussed in the 
previous subsection. 

2.2 Dynamic subgrid-scale modeling for the G-equation 
The main idea of the G-equation is to model flame structure as asymptotically 
thin front. This eliminates the highly nonlinear reaction terms and facilitates model- 
ing for large-eddy simulation. In high Reynolds-number flows, a turbulent premixed 
flame can be viewed as a wrinkled flame “brush” propagating with velocity st • Sev- 
eral previous studies have attempted to derive explicit expressions for st(u / ) (Clavin 
& Williams 1979, Yakhot 1988a, Kerstein and Ashurst 1992). If u ' represents the 
grid- size averaged quantity, this approach is analogous to the original Smagorinsky’s 
subgrid scale model for Navier-Stokes equation in which the eddy viscosity coeffi- 
cient is given a priori. Unfortunately, the existing theoretical and empirical results 
for ST(u f ) do not agree with one another, so that finding the correct functional form 
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of st(u') remains an open question. Even if the question is resolved, there will still 
be a constant to be determined. 

In this section we suggest a new subgrid-scale model for the G-equation based 
on the dynamic modeling principle developed recently (Germano et ai 1991, Moin 
et al 1991). One of the prerequisite conditions for the application of dynamic 
subgrid-scale modeling is that the equation be scale-invariant so that the subgrid 
quantities can be extrapolated from two adjacent scales. The scale-invariance of 
the G-equation has been discussed in the previous studies (Pocheau 1992, Yakhot 
1988b), and was employed in renormalization group theory to derive an explicit 
formula for sx(u') (Yakhot 1988a). We shall skip detailed discussion of this issue. 
We start from the simplest incompressible form of the G-equation; 

f + £(■*» -“FOI < 15 > 

where, although not essential, si is assumed to be constant. Following previous 
works, we define the “grid” filter Q and the “test” filter Q respectively as 

/(x) = J f(x')Q{x,x')dx!, f(x) = Jf(x')Q(x,x')dx', (16) 

where the width of the test filter, A, is larger than that of the grid filter, A. By 
applying the grid filter to (16), we obtain 

% + a?) 

Here both the subgrid scalar flux UjG — UjG and the filtered modulus term |VG| 
need to be modeled. We proceed with applying the test filter, then (17) becomes 

f + s; ( s ' d ) = - sj (P - (18) 

In (17) and (18), it is the filtered modulus term, |VG| that makes the subgrid 
scale modeling of the G-equation difficult compared to other scalar equations. The 
simplest solution is to eliminate this term by applying the test filter to (17) and 
subtract from (18), yielding 

(IS) - (17) = -JJ- (5g - a,d), (19) 

where all the quantities on RHS can now be calculated directly from the large-eddy 
grid solutions. 

We now need to introduce a model to represent the subgrid-scale quantities of 
the G-equation. To this end, we adopt the viewpoint described at the beginning of 
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the subsection, i.e. that, on the large-eddy scale, the turbulent flame brush can be 
represented as a thick front which propagates at speed st . Equations (17) and (18) 
can then be written as 


f + A (a , e)=s> |vc|, (20) 

w + h M) = 5>lv41 ' (21) 

where st and st respectively represent the speed of the flame brush at A and A 
scales. To relate st with the turbulence intensity u', we choose a linear form 


st/sl ~ 1 + C{u /si). 


( 22 ) 


Even if the linear form is not correct, the error may be adjusted by the constant C 
through the dynamic procedure. 

As in the eddy-viscosity model, we further assume u f ~ A|5|, where \S\ = 
\2SijSij\ 1 ! 2 of the large scale strain rate tensor 


Therefore, st and st can be modeled as 

st 


sl 

It 

sl 


= 1 + C G 
= 1 + Cg 


+ St)- 

(23) 

(M), 

(24) 

V SL ) 


(¥■)■ 

(25) 


Substituting (24) and (25) into (20) and (21) and combining with (19) we obtain 



which we wish to use to determine the constant Cg- This is a version of Germano’s 
identity (Germano et al . 1991) for the G-equation. Unlike Germano’s identity used 
in the Navier-Stokes and other scalar equations, however, here we subtract the entire 
equations (17) and (18) instead of treating the subgrid stress terms only, in order 

to eliminate the modulus term |VG| which is difficult to model. Consequently, the 
resulting identity (26) is a single scalar equation for a single unknown parameter Cg, 
rather than the three equations arising from the models for other scalar equations. 
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As in previous work on dynamic subgrid-scale models, the constant Cg is, in 
general, a function of space and time. Therefore, Cg cannot be taken out of the 
test filter, and (26) is an integral equation. However, if the problem of interest has 
at least one homogeneous direction, then Cg can be assumed to be a function of 
the other coordinates and can be removed from the test filter. For example, in the 
case of premixed flame propagating through a channel that is homogeneous in y- 
and ^-directions, Cg = Cg( x,t) so that 



= jvSi H vei-i|-(.,G-s,e), 

(27) 


which is a simpler algebraic equation. 

The modeling proposed in this study is for the simplest constant-density case. 
However, it is anticipated that the same principle can be extended to incorporate 
variable density consideration. The validity of the model is currently under inves- 
tigation for the incompressible G-equation model in homogeneous turbulence. 


3. Future work 

In this study the G-equation model has been applied to several fundamental 
problems relevant to turbulent premixed combustion in the laminar flamelet regime. 
Furthermore, a preliminary dynamic subgrid-scale model for the G-equation has 
been proposed. These ideas need to be further improved to be applied to practical 
high- Reynolds number premixed combustion systems. 

From the standpoint of computational efficiency, the numerical techniques used in 
the present study appear to have a limited application in the practical turbulent re- 
acting flows, partly due to necessity of resolving the abrupt changes in the dependent 
variables across the flame front. It is anticipated that a more efficient discontinuity- 
capturing numerical scheme will greatly reduce the computational cost, thereby 
allowing more extensive parametric studies of fundamental issues such as turbulent 
flame speed. 

As the next step in the application of the large-eddy simulation to combustion, 
the dynamic subgrid-scale model for the G-equation suggested in this study should 
be validated by the direct numerical simulation of the passive G-equation in a 
turbulent flow. If it proves to be successful, then further study is needed to extend 
the model to account for the effects of thermal expansion and variable flame speed. 
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